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accurate  making  it  very  useful  for  routine  modeling  of  2-D  seismic  data.  Its 
principal  problem  is  the  difficulty  in  tracing  rays  in  certain  complex  geo¬ 
logic  models.  Various  calculations  of  synthetic  seismograms  in  laterally 
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ABSTRACT 

A  variety  of  synthetic  seismogram  methods  for  laterally  heterogeneous 
earth  models  have  been  developed  and  tested.  Explicit  finite  difference 
methods  for  elastic  wave  propagation  generally  display  high  accuracy,  but 
are  limited  in  their  application  due  to  extensive  computer  time  and  storage 
requirements.  Implicit  finite  difference  methods  may  be  able  to  realize 
a  factor  of  eight  improvement  in  computer  time  and  a  factor  of  four  improvement 
in  computer  storage  requirements  as  compared  to  explicit  methods.  The 
implicit  methods  of  synthetic  seismogram  calculation  are  less  well  developed 
and  are  more  difficult  in  the  application  of  absorbing  boundary  conditions. 
Approximate  methods  using  the  scalar  wave  equation  may  also  be  useful  for 
preliminary  modeling  of  seismic  data.  A  ray-theoretical  method  using  ray 
tracing  and  disk  ray  theory  (DRT)  synthetic  seismogram  synthesis  has  been 
developed  and  compared  to  calculations  using  reflectivity  (for  1-D)  and 
finite  difference  (for  2-0)  methods.  The  DRT  methcd  is  fast  and  reasonably 
accurate  making  it  very  useful  for  routine  modeling  of  2-D  seismic  data. 

Its  principal  problem  is  the  difficulty  in  tracing  rays  in  certain  complex 
geologic  models.  Various  calculations  of  synthetic  seismograms  in  laterally 
heterogeneous  models  and  comparison  with  homogeneous  models  indicates  the 
importance  of  including  lateral  heterogeneity,  elastic  wave  propagation, 
proper  source  design  and  absorption  in  seismic  modeling. 


INTRODUCTION 


Me  have  been  reviewing  a  variety  of  synthetic  seismogram  techniques 
for  calculation  in  laterally  inhomogeneous  models  for  a  number  of  years. 

In  this  report,  we  review  progress  on  a  variety  of  techniques  including 
disk  ray  theory  (DRT)  based  on  ray  theoretical  methods  and  finite 
difference  (FD)  calculations  which  are  based  on  direct  numerical 
solutions  to  the  seismic  wave  equation  in  two-dimensional  laterally 
heterogeneous  media.  Both  Implicit  and  explicit  formulas  for  finite 
difference  calculations  have  been  studied  and  we  have  utilized  both  the 
acoustic  (scalar)  and  elastic  wave  equations  for  our  finite  difference 
methods.  A  review  of  our  previous  work  in  the  field  of  synthetic 
seismogram  modeling  studies  is  included  in  Department  of  Geosciences, 
Purdue  University,  Report  No.  ONR-1-82  which  was  produced  under  contract 
with  the  Office  of  Naval  Research,  Contract  No.  N00014-75-C-0972.  In 
this  previous  report,  we  described  synthetic  seismogram  modeling 
studies  Including  one-dimensional  and  preliminary  work  with  two-dimen¬ 
sional  models.  In  this  report,  we  show  a  variety  of  synthetic  seismogram 
calculations  primarily  for  laterally  heterogeneous  models  which  illustrate 
the  Importance  of  lateral  Inhomogeneities  in  seismic  wave  propagation 
and  compare  a  variety  of  approaches  to  synthetic  seismogram  calculations. 
Because  these  calculations  can  be  extremely  time  consuming  and  require 
elaborate  computer  programs  and  extensive  computer  time  and  storage, 
we  have  spent  considerable  time  on  attempting  to  find  more  rapid 
methods  which  may  be  approximate,  but  would  provide  for  preliminary 
modeling  In  two  dimensions  such  that  routine  interpretation  could 
utilize  the  vast  and  approximate  techniques  with  final  confirmation 
and  checking  relying  on  the  more  expensive  and  more  accurate  methods. 

Most  of  our  examples  are  for  wave  propagation  and  earth  models 
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corresponding  to  crustal  seismology  studies.  However,  the  techniques 
that  we  describe  are  basically  scale  independent  and  can  be  used  for 
studies  of  seismic  wave  propagation  in  models  corresponding  to  shallow 
engineering  studies  (with  depth  of  penetration  of  several  tens  of 
meters),  to  crustal  studies  (with  depth  of  penetration  of  several  tens 
of  kilometers),  and  to  whole  earth  models  with  propagation  through  the 
entire  earth.  Thus  far  our  calculation  techniques  are  limited  to  two- 
dimensional  propagation  and  flat  earth  models,  although  it  is  conceivable 
that  three-dimensional  methods  and  spherical  earth  geometry  could  be 
utilized  in  the  future. 

EFFECTS  OF  ACOUSTIC  VERSUS  ELASTIC  WAVE  PROPAGATION 
AND  ANELASTICITY  (Q‘] ) 

One  approach  to  approximate  calculation  of  synthetic  seismograms 
in  laterally  inhomogeneous  media  is  to  utilize  an  acoustic  wave 
equation  approach  as  a  preliminary  modeling  technique.  An  acoustic 
finite  difference  formulation  for  example  Is  considerably  faster  than 
the  equivalent  elastic  problem  because  of  the  relative  simplicity  of 
the  scalar  wave  equation  as  compared  to  the  complete  elastic  wave 
equation  in  two  dimensions.  In  utilizing  the  scalar  wave  equation, 
the  complete  nature  of  seismic  waves  Including  the  presence  of  shea<* 
waves  and  the  possibility  of  P-S  converted  phases  is  not  Included.  In 
addition,  some  of  our  finite  difference  methods  utilize  a  perfectly 
elastic  (infinite  Q)  calculation  in  contrast  to  the  more  realistic 
seismic  wave  propagation  in  which  variable  Q  is  Included.  For  these 
reasons,  we  have  calculated  synthetic  seismograms  for  a  one-dimensional 
model  for  which  the  modified  reflectivity  method  can  be  used  to  ill¬ 
ustrate  the  effects  of  acoustic  wave  progagation  and  the  effects  of 
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including  Q  in  the  model  by  comparison  with  a  calculation  in  which 
high  Q  and  elastic  propagation  in  included.  Figure  1  illustrates  a 
seismic  record  section  calculated  for  an  oceanic  crustal  model  in 
which  elastic  wave  propagation  and  high  Q  values  are  present.  The 
model  is  shown  to  the  left  of  the  record  section.  The  upper  layer  Is 
a  200  m  thick  water  layer  in  which  the  shear  velocity  is  set  to  zero 
in  the  program  and  a  density  of  1.0  and  a  P  wave  velocity  of  1.5  km/s 
is  assumed.  The  source  is  within  the  water  layer  at  a  depth  of  100  m. 

The  source  is  approximately  10  Hz.  Relatively  high  Q  values  are  assumed 
in  the  sedimentary  layer  where  a  strong  velocity  gradient  exists  to  a 
depth  of  about  4  km  and  in  the  oceanic  crust  and  upper  mantle  beneath. 

The  seismograms  which  result  from  this  elastic  and  high  Q  calculation 
using  the  modified  reflectivity  method  show  a  complicated  wave  pro¬ 
pagation  with  a  variety  of  P,  multiple  P,  and  P-S  converted  phases  with 
a  range  of  apparent  velocity  from  the  upper  mantle  (8  km/s)  arrivals  to 
a  very  low  apparent  velocity  and  low  frequency  arrival  near  the  end  of 
the  seismograms.  The  high  amplitude/high  frequency  arrival  prominent 
on  all  of  the  seismograms  is  a  water  wave  arrival  with  an  apparent 
velocity  of  approximately  1.5  km/s.  The  low  frequency  arrival  which 
arrives  at  slightly  later  times  than  the  water  wave  is  interpreted  to 
be  a  shear  wave  or  surface  wave  propagating  in  the  sedimentary  layer 
and  coupled  to  water  as  an  acoustic  wave  resulting  in  the  low  apparent 
velocity,  low  frequency,  and  obvious  dispersion.  In  order  to  compare 
the  effects  of  elasticity  versus  acoustic  wave  propagation,  the  identical 
model  as  shown  In  Figure  1  was  utilized  to  calculate  synthetic  seismograms 
assuming  purely  acoustic  wave  propagation.  The  record  section  and  the 
model  are  Illustrated  in  Figure  2.  Acoustic  wave  propagation  in  the 
modified  reflectivity  method  is  provided  for  by  the  transformation 
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suggested  by  Kind  (1976).  The  acoustic  seismograms  shown  in  Figure  2 
are  very  different  from  the  synthetics  for  the  elastic  model.  A  water 
wave  arrival  is  again  present  and  is  in  fact  stronger  than  the  arrivals 
shown  for  the  elastic  model.  However,  most  of  earlier  and  higher 
apparent  velocity  phases  are  either  much  smaller  or  absent.  In  addition, 
the  low  apparent  velocity/low  frequency  arrival  is  also  absent.  Finally, 
increased  numerical  noise  due  to  trapping  of  a  large  percentage  of  the 
energy  in  the  low  velocity  portions  of  the  model  resulting  in  large 
travel  time  arrivals  and  consequent  time  domain  aliasing  in  the  modified 
reflectiy  method  increases  the  numerical  noise  shown  on  the  record 
section.  The  effects  of  Q  are  illustrated  in  Figure  3  in  which  the 
same  elastic  seismic  model  as  shown  in  Figure  1  was  used  to  calculate 
synthetic  seismograms.  However,  the  model  shown  in  Figure  3  includes 
low  Q  values  in  the  sedimentary  layer  and  uppermost  oceanic  crustal 
layer.  For  this  record  section,  all  of  the  same  arrivals  that  were 
included  in  the  record  section  shown  in  Figure  1  are  present,  including 
the  earlier  high  apparent  arrivals,  the  water  wave,  and  the  low  frequency/ 
low  apparent  velocity  arrival.  However,  all  of  these  phases  are  signi¬ 
ficantly  attenuated  and  their  amplitudes  decrease  fairly  rapidly  with 
distance  due  to  the  low  Q  values  in  the  upper  layers. 

These  examples  serve  to  illustrate  the  importance  of  elastic  wave 
propagation  and  the  effects  of  Q  in  synthetic  seismogram  modeling.  The 
elastic  versus  acoustic  assumption  appears  to  be  particularly  important. 
However,  perfectly  elastic  and  acoustic  models  may  be  useful  for 
preliminary  modeling  as  long  as  primary  phases  are  of  principle  import¬ 
ance.  Final  interpretation  of  seismic  data  for  real  earth  structures 
should  include  elastic  wave  propagation  phenomena  and  the  presence  of 
Q  structure. 


EXPLICIT  FINITE  DIFFERENCE  METHODS 


Explicit  finite  difference  methods  involve  a  direct  numerical 
solution  of  the  elastic  wave  equation  in  two  dimensions  for  laterally 
heterogeneous  media.  In  the  explicit  method,  finite  differences  are 
used  to  approximate  spatial  and  temporal  derivatives  in  the  elastic 
wave  equation.  Accuracy  of  the  finite  difference  derivatives  require 
that  approximately  10  points  per  minimum  wavelength  of  the  waves 
propagated  in  the  two-dimensional  model  be  utilized.  Because  the 
displacements  for  time  step  T  +  1  are  calculated  from  two  previous 
time  steps  (T  and  T-l),  the  time  step  must  be  chosen  to  be  very  small 
in  order  to  maintain  stability  of  the  solution  in  this  explicit  method. 
Finite  difference  methods  have  been  presented  by  Boore  (1972),  Alterman 
and  Karal  (1968)  and  Kelly  et  a],.  (1976).  The  calculations  that  we 
illustrate  here  are  from  the  program  developed  by  Mazzella  (1979). 

A  flow  chart  Illustrating  the  explicit  finite  difference  synthetic 
seismogram  calculation  for  acoustic  or  elastic  wave  propagation  is 
shown  in  Figure  4.  An  arbitrarily  complicated  two-dimensional  velocity 
model  is  described  by  a  grid  of  seismic  velocities.  Source  location, 
wavelet  shape  and  peak  frequency  are  designed  which  provide  initial 
displacements  for  the  iterative  finite  difference  calculation.  Displace¬ 
ments  must  be  calculated  for  each  grid  point  location  for  each  time  step 
of  the  model  resulting  in  large  grid  point  and  time  step  loops  in  the 
computer  program.  Figure  5  illunstrates  the  calculation  of  a  displace¬ 
ment  at  a  particular  grid  locatio*  '•Hliz*'  t  the  nine  grid  points  surround 
ing  this  location  at  times  T  and  T-».  The  spatial  derivatives  are  ob¬ 
tained  from  the  displacements  at  the  nine  grid  locations  at  time  equal 
T.  The  time  derivative  Is  obtained  from  the  displacements  at  the 
grid  location  for  times  T  and  T-l.  This  process  must  be  repeated  for 
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all  grid  locations  in  the  model.  The  finite  difference  model  is 
illustrated  in  Figure  6.  The  two-dimensional  velocity  structure  can 
be  arbitrarly  complicated  with  velocities  given  at  each  grid  point. 
Absorbing  boundary  conditions  (Clayton  and  Enquist,  1977,  1980  and 
Enquist  and  Majda,  1977)  help  to  attenuate  the  fictitious  reflections 
obtained  from  artifical  boundaries  of  the  model.  A  free  surface 
boundary  condition  (Ilan,  1975)  is  included  to  account  for  the  free 
surface  conversion.  The  two-dimensional  equations  of  motion  for 
displacement  in  heterogeneous  isotropic  media  and  finite  difference 
approximations  for  the  explicit  case  for  spatial  derivatives  and  time 
derivatives  are  illustrated  in  Figure  7.  The  finite  difference 
approximations  are  straight  forward,  but  are  reasonably  complicated 
and  lead  to  a  long  series  of  steps  in  a  computer  program.  Some  of  the 
practical  difficulties  in  applying  explicit  finite  difference  synthetic 
seismogram  calculations  are  illustrated  in  Figure  8.  Problems  include 
approximations  such  as  two-dimensionality  and  requirements  of  the  grid 
spacing  and  time  step  which  are  related  to  the  maximum  frequency  of 
waves  which  are  to  be  propagated  in  the  model.  These  requirements 
generally  lead  to  large  computer  time  and  storage  needs  for  calculation 
of  realistic  models. 

In  order  to  test  our  finite  difference  program,  we  have  calculated 
synthetic  seismograms  using  the  modified  reflectivity  method  and  the 
finite  difference  program  for  a  layer  over  a  half-space  model.  The 
resulting  record  section  (Figure  9)  for  vertical  component  seismograms 
show  good  comparison  for  primary  P  wave  arrivals,  multiples,  and  P-S 
conversions.  Some  discrepancy  in  surface  wave  arr.vals  is  apparent 
due  to  the  phase  velocity  range  Included  in  the  modified  reflectivity 
method  calculation  and  due  to  the  presence  of  some  numerical  dispersion 
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in  the  finite  difference  calculation  for  the  relatively  short  wave¬ 
length  surface  waves  due  to  their  low  propagation  velocity.  A  more 
detailed  comparison  of  the  seismograms  for  the  finite  difference  and 
modified  reflectivity  methods  is  shown  in  Figure  10  for  distances 
of  10,  20,  30  and  40  km.  The  modified  reflectivity  seismograms  for 
these  distances  are  illustrated  immediately  beneath  the  finite  difference 
seismograms.  Good  amplitude  and  waveform  comparison  is  seen  for  all 
phases  except  the  surface  wave. 

The  finite  difference  elastic  wave  equation  method  represents  an 
accurate  numerical  technique  for  the  calculation  of  two-dimensional 
synthetic  seismograms  in  laterally  heterogeneous  media.  At  present, 
computer  limitations  preclude  its  use  in  routine  modeling  of  seismic 
data.  For  most  models  of  geophysical  and  geological  interest,  the 
finite  difference  calculation  of  realistic  synthetics  with  the  finite 
difference  program  would  require  computer  time  on  the  order  of  several 
hours  using  the  fastest  computers  available  (Cray  1  and  Cyber  205) 
and  in  addition  would  require  core  storage  of  greater  than  a  million 
words.  Therefore,  we  consider  the  finite  difference  approach  at  the 
present  time  to  be  primarily  of  interest  for  model  studies  in  which 
various  kinds  of  wave  propagation  phenomena  can  be  studied  with 
seismograms  for  characteristic  models.  In  addition,  the  finite 
difference  program  synthetics  represent  an  important  check  on  faster, 
but  more  approximate  techniques. 

DISK  RAY  THEORY  SYNTHETIC  SEISMOGRAM  METHOD 

A  disk  ray  theory  (DRT)  synthetic  seismogram  method  for  approximate 
calculation  of  seismic  wave  propagation  in  laterally  heterogeneous  elastic 
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and  anelastic  media  has  been  developed  and  tested.  The  disk  ray  theory 
method  follows  the  theory  developed  by  Wiggins  and  Madrid  (1974)  and 
Wiggins  (1976).  The  technique  utilizes  ray  tracing  and  complex  two- 
dimensional  velocity  and  Q  models  to  provide  the  necessary  time  distance 
and  amplitude  data  to  construct  synthetic  seismograms  by  a  summation 
of  the  contribuitons  of  the  various  raypaths.  A  description  of  the 
technique  and  some  examples  of  its  application  are  shown  in  a  paper 
entitled  "An  Example  of  Application  of  Synthetic  Seismogram  Modeling 
in  Two  Dimensions"  by  C.S.  Chiang  and  L.W.  Braile.  A  copy  of  this 
paper  is  included  in  the  Appendix  to  this  report. 

The  principal  advantage  of  the  DRT  synthetic  seismogram  method  is 
speed.  It  is  relatively  accurate  and  synthetics  can  be  quickly  constructed 
for  two-(or  even  three-)  dimensional  models  for  which  ray  tracing  can 
be  performed.  The  method  can  include  primary  P  and  S  phases  as  well 
as  multiples  (Figure  11).  Its  principal  limitation  is  that  guided  or 
interference  propagation  type  of  phases  cannot  be  included.  Surface 
waves  are  not  accounted  for  and  some  geologic  models  are  difficult  or 
impossible  to  trace  rays  through  so  that  the  necessary  ray  trace  input 
may  not  be  available  for  certain  physically  important  phases.  An 
example  of  a  DRT  calculation  is  illustrated  for  the  model  shown  in 
Figure  12.  The  seismic  raypaths  are  traced  using  a  ray  tracing  program 
for  the  phases  shown.  Amplitude,  travel  time,  raypath  parameter  and 
distance  date  are  then  used  to  construct  synthetic  seismograms  by  the 
DRT  method.  For  this  relatively  simple,  but  important  laterally 
inhomogeneous  fault  model,  the  seismic  rays  can  be  traced  with  high 
accuracy  and  efficiency  to  produce  the  data  for  the  DRT  seismogram 
synthesis.  This  model,  however,  illustrates  one  of  the  problems  with 
the  DRT  method.  The  seismic  phases  shown  in  Figure  12  corresponding 
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to  the  propagation  at  the  bottom  of  the  prominent  fault  (Phases 
Prefl-2  and  Ph-2^  cannot  be  traced  directly  in  this  model.  If  one 
attempts  to  trace  these  raypaths  from  the  source  to  the  bottom  of  the 
fault  plane  (corner)  the  raypaths  do  not  refract  by  Snell's  Law  to  the 
received  locations  as  shown.  In  fact,  the  lower  corner  of  the  fault 
represents  a  diffractor.  However,  we  know  from  wave  theoretical 
principals  and  from  a  finite  difference  calculation  that  the  waves 
corresponding  to  these  raypaths  are  physically  important  and  are  present 
in  the  seismic  wave  propagation  in  this  model.  Therefore,  we  were  able 
to  generate  these  phases  by  a  two  step  ray  tracing  procedure.  All 
other  phases  were  ray  traced  in  the  normal  method  from  source  to  reflecting 
or  refracting  point  and  then  to  receiver.  However,  for  these  two  phases, 
rays  were  traced  from  the  source  to  the  diffraction  point  (corner  at 
the  bottom  of  the  fault).  Then,  an  artifical  source  was  located  at 
this  corner  and  caused  to  generate  the  phases  prefi_2  and  Ph_2-  The 
travel  time  and  amplitude  factors  for  this  artifical  diffraction  source 
were  then  combined  with  the  raypaths  from  the  source  to  the  corner  to 
obtain  the  time  distance  and  amplitude  factors  necessary  for  DRT 
synthesis  for  these  complicated  raypaths.  In  this  way,  the  correct 
wave  propagation  phenomena  were  simulated.  However,  the  amplitude 
factors  which  resulted  from  this  calculation  were  too  large  of 
amplitude  for  these  two  phases  from  the  point  diffractor.  By  assuming 
that  the  effective  reflection  coefficient  of  the  point  diffractor  at 
the  corner  of  the  fault  was  approximately  0.3,  seismograms  from  the 
DRT  method  which  compare  favorably  with  finite  difference  synthetic 
were  produced  (Figure  13).  The  finite  difference  and  DRT  synthetics 
for  the  fault  model  (INFL-U)  shows  excellent  comparison  for  all  of 
the  primary  arrivals.  Later  arrivals  (multiple  P  waves,  P-S  conversions 


and  surface  waves)  included  in  the  finite  difference  calculation  are 
not  included  in  the  DRT  synthesis.  In  addition,  the  wide  angle  re¬ 
flection  from  the  upper  interface  (Phase  Pref|_-|)  has  too  large  of 
an  amplitude  at  large  distances  for  the  DRT  method.  This  is  a  result 
of  the  limitation  of  ray  theory  for  grazing  incidence  wave  propagation 
at  low  frequencies  for  thin  layers.  However,  this  is  not  severe 
limitation  of  the  DRT  method  because  such  arrivals  are  rarely  of 
significance  anyway  because  of  the  lower  apparent  velocity  of  such 
phases  and  the  presence  of  attenuation  which  normally  decreases  their 
amplitudes  fairly  rapidly.  A  more  detailed  comparison  of  the  finite 
difference  and  disk  ray  theory  synthetics  for  the  fault  model  (INFL-11) 
is  shown  in  Figure  14.  Amplitude  and  wave  form  character  of  phases 
of  interest  are  nearly  identical  for  the  FD  and  DRT  method.  Therefore, 
we  conclude  that  the  DRT  method  provides  a  very  useful  approach  to 
preliminary  modeling  of  seismic  data  in  two-dimensional  or  even  three- 
dimensional  media.  The  technique  is  fast,  reasonably  accurate  and  is 
limited  only  due  to  the  difficulty  of  tracing  rays  for  certain  phases 
and  certain  geological  models.  We  suggest  that  this  technique  would 
provide  a  very  useful  modeling  method  in  which  routine  ray-tracing 
and  DRT  synthesis  can  provide  comparison  with  observed  seismic  data 
for  two-dimensional  velocity  structures.  For  final  confirmation  or 
where  the  technique  is  not  appropriate,  finite  difference  or  other 
wave  theortical  techniques  will  still  be  necessary. 

FINITE  DIFFERENCE  SYNTHETIC  SEISMOGRAM  CALCULATIONS 
FOR  DISLOCATION  SOURCE  MODELS. 

Explicit  finite  difference  synthetic  seismogram  calculations 
program  for  dislocations  source  models  has  been  presented  by  Espindola 
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(1979).  The  technique  follows  similar  methods  as  illustrated  previously 
in  this  report  and  are  also  similar  to  those  presented  by  Mazzella  (1979). 
However,  in  this  method  a  skew  symmetric  coordinate  system  (Figure  15) 
is  utilized  to  provide  for  a  distributed  dislocation  source  which 
simulates  an  earthquake.  A  dislocation  with  arbitrary  slip  function  and 
source  time  history  is  caused  to  occur  at  grid  points  corresponding  to 
the  source  fault  plane  location.  Finite  difference  calculations  allow 
for  propagation  of  seismic  waves  away  from  this  distributed  source. 

The  skew  symmetric  coordinate  system  is  required  in  order  to  provide  for 
parallel  1  ism  of  the  grid  spacing  with  the  source  fault  plane.  This 
program  therefore,  has  an  additional  advantage  over  the  previously 
described  finite  difference  program  in  that  a  variety  of  more  realistic 
source  phenomena  can  be  studies  using  the  finite  difference  method. 

The  program  of  the  method  still  utilizes  the  explicit  approach  and 
is  limited,  therefore  requiring  certain  grid  point  and  time  step 
conditions,  necessitating  large  computer  time  and  storage  and  is  also 
limited  to  a  two-dimensional  approximation  such  that  the  source  fault 
plane  is  assumed  to  extend  in  the  direction  perpendicular  to  the  XZ 
axis.  Examples  of  synthetic  seismogram  calculations  using  this  dis¬ 
location  finite  difference  method  are  shown  for  a  homogeneous  half 
space  model  (Figure  16)  and  a  complicated  laterally  heterogeneous 
velocity  structure  model  (Figure  17).  A  variety  of  source  slip 
functions  and  source  time  histories  for  the  identical  fault  geometry 
(Figure  18)  are  also  Illustrated.  Figure  19  shows  the  seismic  ray 
paths  which  are  produced  by  the  distributed  source  dislocation  for 
a  simple  displacement  time  history.  The  prominent  effects  of  two- 
dimensionality  can  be  illustrated  by  comparing  seismograms  calculated 
for  the  homogeneous  and  the  laterally  heterogeneous  models.  Figure  20 
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shows  vertical  and  radial  component  seismograms  (WWSSN  response) 
calculated  for  model  S-2  with  a  ramp  function  displacement  time  history. 
Comparison  of  selected  seismograms  for  the  homogeneous  and  laterally 
heterogeneous  calculations  (models  S-1A  and  S-2A)  are  shown  in  Figure 
21.  Both  short  period  and  long  period  WWSSN  seismograms  responses 
are  shown.  The  effects  of  inhomogeneity  are  prominent  primarily  in 
later  arriving  phases  on  Z  and  R  components  of  model  S-2A  and  in  the 
large  amplitude  phase  which  is  present  primarily  on  the  radial  com¬ 
ponent  for  model  S-2A.  Long  period  records  also  show  considerable 
difference  in  wave  form  due  to  the  wave  propagation  and  the  complex 
structure.  The  Identical  source  is  used  for  both  of  these  models. 

Comparison  of  the  effects  of  differing  source  functions  is 
illustrated  in  the  seismograms  shown  in  Figure  22.  All  of  these 
synthetic  seismograms  are  calculated  for  the  identical  homogeneous 
half  space  model  and  Identical  fault  geometry  (Figure  16  and  Figure 
18).  However,  the  source  time  history  and  slip  functions  illustrated 
in  Figure  18  were  used  for  calculations  S-1A,  S-1B  and  S-1C.  The 
ramp  source  time  history  (S-1A)  produces  considerably  more  high 
frequency  wave  propagation  and  more  complex  phase  arrivals  than  for 
the  smoother  slip  function  and  time  history  calculations. 

These  calculations  serve  to  illustrate  the  utility  of  the  finite 
difference  method  for  application  to  earthquake  problems  in  which  the 
source  function  is  of  significance  and  also  to  demonstrate  the  Importance 
of  lateral  heterogeneity  in  wave  propagation  for  models  of  complex 
geologic  structure. 
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IMPLICIT  FINITE  DIFFERENCE  SYNTHETIC  SEISMOGRAMS  CALCULATIONS 

We  have  pursued  the  use  of  implicit  numerical  techniques  in  finite 
difference  formulations,  in  hopes  of  reducing  both  computer  time  and 
storage  requirements  of  existing  predominately  explicit  formulations. 
Implicit  techniques,  which  involve  centered-differences  calculations 
at  each  time  step,  have  two  potential  advantages  over  explicit  formulations. 
Firstly,  some  of  the  implicit  formulations  are  numerically  stable  for  any 
step  size  and  grid  spacing.  This  characteristic  allows  the  length  of 
the  time  steps  to  be  increased,  thus  decreasing  the  actual  number  of 
time  steps  without  any  regard  for  stability.  One  must  use  caution  in 
increasing  the  size  of  the  time  step,  however,  to  insure  that  the 
accuracy  is  maintained.  Secondly,  another  potential  advantage  of 
implicit  formulations  Is  that  some  of  the  techniques  are  more  accurate 
than  their  explicit  couterparts.  Thus,  a  coarser  grid  spacing  along 
with  a  corresponding  increase  in  time  step  size  is  afforded  while 
maintaining  accuracy  and  stability  for  certain  implicit  formulations 
with  higher  accuracy.  In  general,  our  experimentation  indicates  a 
trade  off  between  stability  and  accuracy  of  implicit  formulations. 

If  both  time  steps  and  grid  spaclngs  can  be  Increased,  significant 
savings  in  computer  time  and  storage  requirements  may  be  possible 
with  implicit  methods.  For  the  most  part  our  implicit  studies  have 
utilized  the  simpler  scalar  (acoustic)  wave  equations.  Experience 
with  these  techniques  indicates  that  further  application  to  the  elastic 
wave  equation  may  be  warranted.  However,  additional  problems,  such  as 
the  difficulty  of  applying  absorbing  boundary  conditions  to  the  more 
complicated  implicit  formulas  for  the  elastic  case  may  diminish  some  of 
the  anticipated  benefits  of  implicit  formulas. 


The  principle  disadvantage  of  implicit  finite  difference  formulations 
are  primarily  concerned  with  their  complexity.  Implicit  calculations  are 
not  as  straight  forward  as  their  explicit  couterparts  (Figure  23  and 
Figure  24).  Implicit  difference  equations  for  two  or  three-dimensional 
problems  have  to  be  broken  down  corresponding  to  two  or  three  elements 
using  either  an  alternating  direction  implicit  (ADI)  or  a  locally  one 
dimensional  (LOD)  spliting  technique.  Each  of  these  parts  must  be 
solved  inversely  using  a  tri -diagonal  matrix  equation  solver.  These 
equations  are  particularly  difficult  to  formulate  for  a  fully  absorbing 
boundary  condition  at  the  edges  and  bottom  of  the  model,  and  for  a  free 
surface  boundary  condition  at  the  top  of  the  model. 

We  have  compared  both  with  accuracy  calculations  and  actual  implicit 
finite  difference  programs,  a  variety  of  implicit  finite  difference 
techniques  (Figure  25).  The  different  formulas  illustrated  in  Figure  25 
have  different  accuracy  and  stability  conditions  and  the  complexity  of 
the  formulas  leads  to  differing  degrees  of  difficulty  in  applying 
boundary  conditions  source  functions  and  other  practical  aspects  of 
producing  a  workable  finite  difference  program  for  two-dimensional 
seismic  wave  propagation  problems. 

One  of  our  most  significant  experimental  results  has  come  from 
comparing  four  different  implicit  formulations  of  the  scalar  wave 
equation  with  an  explicit  formulation  for  accuracy.  Our  accuracy 
criterian  is  derived  from  the  numerical  dispersion  studies  published 
by  Alford  et  al .  (1974)  and  Emerman  et  al .  (1982).  However,  we  have 
investigated  the  numerical  dispersion  effects  of  all  the  formulas  shown 
in  Figure  25.  Some  of  the  results  are  illustrated  in  Figures  26  through 
29.  While  some  of  Implicit  formulations  prove  to  be  less  accurate  than 
the  explicit  case,  the  equations  proposed  by  Fairweather  and  Mitchell 
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prove  to  be  far  superior.  For  example,  one  can  see  from  the  explicit 
accuracy  diagram  shown  in  Figure  26  that  in  order  to  maintain  reasonable 
accuracy  in  an  explicit  solution,  one  must  use  approximately  ten  points 
per  wave  length  (1  *  0.1)  in  order  to  prevent  significant  numerical 

la 

dispersion.  This  ten  points  per  wave  length  requirement,  in  addition 

to  the  small  times  step  required  for  stability,  makes  the  explicit 

case  extremely  costly  in  computer  time  and  storage.  The  implicit 

finite  difference  formula  presented  by  Emerman,  Schmidt  and  Stephen 

(1982)  does  not  have  the  stability  requirement,  however,  one  can  see 

from  Figure  27  that  a  small  time  step  is  again  required  in  order  to 

maintain  reasonable  accuracy.  Further,  the  accuracy  of  the  solution 

decreases  with  increasing  time  step  (increasing  p  value)  making  this 

formulation  rather  undesirable  because  although  the  time  step  can  be 

increased  and  stability  still  maintained  thus  improving  on  computer 

time,  the  grid  spacing  must  be  kept  extremely  small  (greater  than  10 

points  per  wave  length)  in  order  to  maintain  accuracy.  However,  the 

numerical  dispersion  accuracy  curves  (Figure  28  and  Figure  29)  for  the 

Fairweather  and  Mitchell  (1965)  formula  show  that  this  formula  is  far 

superior.  For  this  case,  minimal  dispersion  is  observed  for  as  few  as 

four  to  five  grid  points  per  wave  length  (]_  =  0.2  to  0.25)  and  for  a 

G 

variety  of  time  step  conditions  (p  =  0.1  to  p  *  0.7),  thus  one  can 
maintain  reasonable  accuracy  with  the  Fairweather  and  Mitchell  formula 
with  approximately  twice  as  coarse  of  a  grid  spacing  as  for  the  explicit 
or  other  implicit  methods.  For  a  two-dimensional  model  these  conditions 
result  In  a  factor  of  four  savings  in  computer  storage  and  of  eight 
savings  in  computer  time,  because  the  time  step  can  also  be  increased  as 
the  grid  spacing  is  increased.  The  Fairweather  and  Mitchell  formula, 
however,  does  have  a  stability  limit,  such  that  the  time  step  can  only 
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be  increased  to  p  *  0.7.  Nevertheless,  this  potential  savings  in 
computer  time  of  approximately  a  factor  of  eight  and  a  savings  of 
computer  storage  requirements  of  approximately  a  factor  of  four  make 
the  Falrweather  and  Mitchell  implicit  formula  worth  investigating 
further. 

We  have  developed  working  programs  using  the  implicit  finite 
difference  formulation  for  the  acoustic  (scalar)  wave  equation.  We 
have  used  both  the  Emerman,  Schmidt  and  Stephen  implicit  formulas  as 
well  as  the  Fairweather  and  Mitchell  formulas.  Although  further 
analysis  is  required  and  detailed  comparisons  of  accuracy,  stability 
and  calculation  time  between  the  two  implicit  methods  and  the  explicit 
approach  are  required,  we  believe  that  the  Fairweather  and  Mitchell 
implicit  finite  difference  method  will  be  useful  for  further  application 
to  synthetic  seismogram  calculation  in  two-dimensional  models.  We  have 
succeeded  in  applying  relatively  accurate  non-absorbing  boundary 
conditions,  at  the  edges  of  implicit  models  (Figure  30).  However,  the 
higher  order  absorbing  boundary  conditions  need  to  be  applied  to  the 
Fairweather  and  Mitchell  formulas  and  these  implicit  techniques  need  to 
be  applied  to  the  elastic  wave  equation  case  in  order  to  fully  realize 
the  benefits  of  implicit  difference  methods  in  two-dimensional  synthetic 
siesmogram  calculation. 
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FIGURE  CAPTIONS 


Figure  1.  Synthetic  seismograms  calculated  for  a  plane  layered  shallow 

oceanic  model.  The  calculation  includes  the  effects  of  elastic 
wave  propagation  in  the  oceanic  crustal  layers  beneath  a  200 
meter  depth  and  acoustic  wave  propagation  in  the  upper  water 
layer  (200  meters  thick).  The  frequency  of  the  source  is 
approximately  10  Hz  and  the  source  depth  is  100  meters.  High 
Q  values  (Qp  =  1,000,  500,  1,000  are  assumed  for  the  water 

layer,  sedimentary  layer  and  oceanic  crustal  and  upper  mantle 
layers,  respectively).  Amplitudes  are  scaled  for  plotting  by 
multiplying  the  amplitude  of  each  seismogram  by  the  distance. 

Figure  2.  Synthetic  seismograms  calculated  for  an  acoustic  layered  half- 
space.  The  layer  parameters  are  the  same  as  for  Figure  1. 
Acoustic  wave  propagation  is  simulated  in  both  the  water 
layer  and  the  oceanic  crustal  and  upper  mantle  layers.  High 
Q  values,  as  in  Figure  1,  are  assumed.  The  acoustic  propagation 
in  the  reflectivity  program  is  simulated  by  the  transformation 
described  by  Kind,  1976.  Some  numerical  noise  in  the  calculation 
is  evident  throughout  the  record  section.  The  amplitudes  for 
this  seismogram  are  scaled  to  one-half  of  the  amplitudes  shown 
on  Figure  1 . 

Figure  3.  Synthetic  seismograms  calculated  for  the  identical  oceanic 
model  as  in  Figure  1,  except  that  low  Q  values  are  included 
in  the  sedimentary  layer  and  the  upper  oceanic  crystal  layer. 

An  elastic  half-space  beneath  the  acoustic  water  layer  (thick¬ 
ness  of  200  meters)  is  assumed.  Amplitude  scaling  is  identical 
to  that  shown  on  Figure  1. 

Figure  4.  Flow  chart  for  the  calculation  of  finite  difference  synthetic 
seismograms  for  acoustic  or  elastic  earth  models  using  the 
explicit  finite  difference  formulation. 

Figure  5.  Illunstration  of  the  computational  technique  for  spatial  and 
temporal  finite  difference  derivatives  at  a  time  T+l  using 
the  surrounding  displacement  grid  points  at  times  T  and  T-l. 

Figure  6.  Schematic  diagram  illustrating  the  calculation  of  finite 

difference  synthetic  seismograms  for  a  two-dimensional  earth 
model.  The  velocity  structure  is  approximated  by  a  grid  in 
which  the  seismic  velocity  is  given  at  each  grid  location. 

Seismic  waves  are  caused  to  propagate  by  a  finite  difference 
time  step  calculation  away  from  an  arbitrarily  located  source 
point.  A  free  surface  boundary  condition  and  absorbing 
boundaries  at  the  fictitious  edges  of  the  model  are  included. 
Seismometer  locations  record  the  displacement  time  history  for 
the  wave  propagation  in  the  two-dimensional  model. 


Figure  7.  Two-dimensional  equations  of  motion  for  displacement  in 

heterogeneous  isotropic  media  and  finite  difference  approximations 
for  the  explicit  finite  difference  formulation  of  second  order 
spatial  derivatives,  cross-product  spatial  derivatives  and  second 
order  time  derivatives  for  approximations  to  the  two-dimensional 
equation  of  motion. 

Figure  8.  List  of  some  of  the  practical  aspects  of  explicit  finite 
difference  synthetic  seismogram  calculation. 

Figure  9.  An  example  of  finite  difference  synthetic  seismogram  calculation 
in  a  laterally  homogeneous  media.  Upper  record  section  shows 
synthetic  seismograms  calculated  for  a  layer  over  a  half-space 
model  using  the  modified  reflectivity  method.  Lower  record 
Section  shows  synthetic  seismograms  calculated  by  the  finite 
difference  technique  for  the  identical  model.  The  model 
consisted  of  a  2  km  thick  4.5  km/s  layer  overlying  a  6  km/s 
half-space.  In  the  upper  diagram,  amplitudes  are  scaled  by 
distance  squared.  In  the  lower  record  section,  amplitudes 
are  scaled  by  distance  to  the  one  power.  The  difference  in 
amplitude  scaling  corresponds  to  the  difference  in  assumptions 
of  the  two  methods.  In  the  modified  reflectivity  method, 
a  point  source  in  a  one-dimensional  model  represents  three- 
dimensional  spreading  of  the  wave  energy.  In  the  finite 
difference  model,  however,  a  two-dimensional  approximation 
is  assumed  and  the  source  is  actually  a  line  source. 

Some  differences  in  the  later  portions  of  the  record 
section  are  observable  due  to  the  fact  that  the  surface 
waves  were  not  exactly  included  in  this  particular  calculation 
using  the  modified  reflectivity  method  and  also  due  to  the 
fact  that  because  of  the  choice  of  grid  spacing  and  time 
step  parameters  in  the  finite  difference  calculation,  some 
numerical  dispersion  is  included  for  shear  wave  and  surface 
wave  propagation.  However,  the  compress ional  waves  for 
primary  and  multiple  reflections  are  correctly  included  in 
both  methods  and  a  comparison  of  the  two  synthetic  seismogram 
record  sections  indicates  that  the  finite  difference  seismograms 
are  nearly  identical  to  the  modified  reflectivity  method 
sei smograms . 

Figure  10.  Comparison  of  finite  difference  and  modified  reflectivity 
synthetic  seismograms  for  the  layer  over  a  half-space  model 
described  in  Figure  9.  These  seismograms  illustrate  a  detailed 
comparison  of  the  synthetics  at  distances  of  10,  20,  30  and 
40  km  from  the  source. 
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Figure  11.  Synthetic  seismograms  calculated  for  the  modified  reflectivity 
and  disk  ray  theory  methods  for  a  layer  over  a  half-space 
model  including  primary  congressional  waves,  first  order 
multiples  and  selected  P-S  conversions.  The  DRT  synthetics 
are  shifted  slightly  to  the  left  of  their  proper  distance 
for  convenient  display  and  comparison  with  the  MRM  seismograms. 
Amplitudes  of  the  seismograms  are  scaled  by  distance  to  the 
one  power  for  convenient  plotting. 

Figure  12.  Seismic  raypaths  for  DRT  calculation  for  fault  model  INFL-11. 

Figure  13.  Seismic  record  sections  for  finite  difference  and  disk  ray 

theory  synthetic  seismograms  for  fault  model  INFL-11.  Seismic 
phases  are  identified  for  the  raypaths  shown  in  Figure  12. 

Finite  difference  synthetic  seismograms  are  scaled  by  multiplying 
by  distance  to  the  one  power,  whereas  the  disk  ray  theory 
synthetics  are  scaled  by  multiplying  by  the  distance  squared 
to  account  for  the  difference  in  line  source  for  the  finite 
difference  program  versus  point  source  for  the  disk  ray 
theory  program.  Multiples  and  P-S  converted  phases  are  not 
included  for  the  disk  ray  theory  calculation,  but  are  included 
in  the  finite  difference  model. 

Figure  14.  Comparison  of  finite  difference  and  disk  ray  theory  synthetics 
for  fault  model  INFL-11  for  distances  10,  20,  30,  40  and 
50  km.  Phase  notation  is  identical  to  that  shown  on  the 
raypath  diagram  (Figure  12). 

Figure  15.  Schematic  diagram  illustrating  the  grid  geometry  for  finite 
difference  model  calculations  in  a  skew  coordinate  system 
with  a  dislocation  (earthquake)  source. 

Figure  16.  Homogeneous  half-space  model  for  skew  coordinate  system 
finite  difference  calculation.  The  source  is  a  reversed 
fault  beneath  seismometer  location  number  6.  The  source 
time  history  is  shown  as  an  inset  in  the  model.  A  unilateral 
source  is  assumed  with  fracturing  beginning  at  the  lower 
end  of  the  fault  plane. 

Figure  17.  A  laterally  inhomogeneous  earth  model,  representing  a 
subduction  zone,  with  identical  geometry  and  source  to 
the  homogeneous  model  shown  in  Figure  16. 

Figure  18.  Diagram  illustrating  the  fault  geometry  slip  functions  and 
source  time  histories  for  a  variety  of  calculations  using 
the  skew  symmetric  coordinate  system  finite  difference 
synthetic  seismogram  program. 
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Figure  19.  Schematic  diagram  illustrating  the  seismic  raypaths  which 
produce  prominent  arrivals  from  a  dipping  fault  source  for 
a  ramp  function  displacement  in  the  homogeneous  half-space 
model  (Figure  16).  Dashed  raypaths  indicate  S  wave  pro¬ 
pagation  and  solid  raypaths  indicate  congressional  wave 
propagation.  Unprimed  phase  notation  indicates  arrivals 
from  the  starting  phase  of  the  source  time  history,  whereas 
primed  phases  are  due  to  the  stopping  phase  two  seconds 
after  the  initiation  of  fracturing.  Subscript  1  indicates 
arrivals  from  the  point  of  initiation  of  fracture  at  the 
bottom  of  the  fault  plane.  Subscript  2  indicates  arrivals 
propagating  from  the  top  end  of  the  fault  plane. 

Figure  20.  Synthetic  seismograms  calculated  for  model  S2  (Figure  17) 
for  vertical  and  radial  component  displacements.  The 
seismograms  have  been  convolved  with  a  seismograph  response 
function  corresponding  to  a  short  period  WWSSN  response. 
Seismograms  are  plotted  as  a  function  of  distance  from  the 
zero  point  directly  above  the  point  of  initiation  of  the 
earthquake  source.  Small  numbers  adjacent  to  the  seismograms 
indicate  the  seismometer  position.  Relative  amplitudes  are 
constant  for  all  seismograms  and  can  be  compared  with  the 
relative  amplitude  scale  shown  at  the  bottom  of  the  record 
section. 

Figure  21.  Comparison  of  short  period  and  long  period  seismograms  for 
models  SI A  and  S2A  at  a  distance  of  -6.1  km  corresponding 
to  seismometer  number  6.  Both  vertical  (SPZ  and  LPZ)  and 
radial  (SPR  and  LPR)  seismograms  are  shown.  The  long  period 
records  are  scaled  with  a  difference  amplitude  scale  as 
indicated  at  the  bottom  of  the  diagram. 

Figure  22.  Comparison  of  SPZ  seismograms  for  the  homogeneous  half-space 
fault  model  (Figure  16)  with  the  three  difference  source 
functions  (1A,  IB,  1C)  illustrated  in  Figure  18.  The  phase 
notation  is  the  same  as  illustrated  in  Figure  19.  Dotted 
lines  indicate  portions  of  seismograms  with  an  enhanced 
amplitude  scale  to  see  small  arrivals.  The  amplitude  scales 
differ  for  each  of  the  record  sections  and  are  scaled  according 
to  the  relative  amplitude  scale  bar  shown  at  the  bottom  of 
each  record  section. 

Figure  23.  Equations  used  in  calculating  acoustic  finite  difference 
synthetic  seismograms  using  the  scalar  wave  equation.  An 
explicit  and  an  implicit  finite  difference  form  are  illustrated 
as  well  as  second  order  finite  difference  operators. 

Figure  24.  Schematic  diagram  illustrating  the  method  of  propagation 

of  pressures  at  time  N+l  from  times  N  and  N-l  using  the  explicit 
finite  difference  form  as  well  as  a  similar  diagram  illustrating 
the  same  calculation  using  the  implicit  form. 
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Figure  25.  Five  different  finite  difference  schemes  for  explicit  and 
implicit  finite  difference  calculations  using  the  scalar 
(acoustic)  wave  equation. 

Figure  26.  Accuracy  diagram  for  the  explicit  acoustic  wave  equation. 

The  normalized  phase  velocity  is  a  measure  of  the  degree 
of  numerical  dispersion  present  in  the  calculation  for  a 
variety  of  grid  spacings.  The  value  of  1/G  of  0.1  which 
corresponds  to  10  points  per  wavelength  is  normally  considered 
appropriate  for  explicit  finite  difference  calculations  to 
maintain  a  negligible  amount  of  numerical  dispersion.  The 
calculation  is  for  an  angle  of  incidence  with  respect  to  the 
grid  axes  of  0°.  The  P  values  are  relative  to  the  time  step 
of  the  calculation  where  P  =  Cq  x  aT/aX. 

Figure  27.  Accuracy  diagram  for  the  implicit  formula  given  by  Emerman, 
Schmidt  and  Stephen  (1982). 

Figure  28.  Accuracy  diagram  for  the  implicit  finite  difference  equation 
given  by  Fairweather  and  Mitchell  (1965). 

Figure  29.  Accuracy  diagram  for  the  implicit  finite  difference  formula 
given  by  Fairweather  and  Mitchell  (1965).  In  this  case, 
the  P  value  is  kept  constant  at  0.7  and  the  various  angles 
of  incidence  from  0  to  45°  are  assumed. 

Figure  30.  Displacement  time  history  diagrams  (snapshots)  of  an  implicit 
finite  difference  calculation  in  a  homogeneous  media  in  two 
dimensions.  The  two-dimensional  space  and  the  displacements 
are  illustrated  in  a  perspective  view.  The  edges  of  the 
perspective  diagram  correspond  to  the  artifical  edges  of  the 
model.  A  non-reflecting  boundary  condition  is  assumed  and 
relatively  small  amounts  of  energy  are  reflected  back  from 
the  boundaries. 
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Figure  24. 
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ABSTRACT 


Synthetic  seismograms  for  laterally  inhomogeneous  velocity  and  Q 
structures  can  be  calculated  by  the  disk  ray  theory  (DRT)  method.  Ray¬ 
tracing  through  the  model  provides  travel -time,  distance,  slowness  and 
ray  amplitude  data  which  are  combined  to  synthesize  seismic  record  sections 
by  DRT.  The  amplitude  factors  include  geometrical  spreading,  attenua¬ 
tion  due  to  Q'1,  reflection  and  transmission  coefficients  and  the  free 
surface  conversion  coefficient.  Comparison  of  DRT  synthetic  seismograms 
with  those  calculated  by  the  reflectivity  method  for  a  layer  over  a  half 
space  model  shows  excellent  agreement  in  both  amplitude  and  waveform 
character  of  all  arrivals.  Direct,  reflected,  P  to  S  converted,  multiply 
reflected  and  head  waves  may  be  synthesized  by  the  DRT  method.  In  an 
application  of  two-dimensional  synthetic  seismogram  modeling  using  the 
DRT  method,  synthetic  seismogram  record  sections  were  calculated  for 
a  complicated  geologic  model  of  the  crust  in  the  eastern  Snake  River 
Plain,  Idaho  and  compared  with  observed  seismic  data.  In  this  application, 
we  find  that  a  small  modification  of  the  Snake  River  Plain  velocity  model 
which  was  interpreted  from  travel -time  data  produces  seismograms  which 
compare  well  with  the  amplitude  and  waveform  character  of  the  observed 
seismograms. 
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INTRODUCTION 


Modern  crustal  seismic  refraction  and  reflection  data  recorded  at 
relatively  small  station  spacings  and  of  improved  quality  result  in  an 
increased  need  for  two-dimensional  synthetic  seismogram  modeling  tech¬ 
niques  for  interpretation.  Both  wave- theoretical  and  ray-theoretical 
methods  have  been  used  for  the  calculation  of  synthetic  seismograms. 

Wave  methods  are  generally  restricted  to  one-dimensional  velocity  models, 
although  finite-element  and  finite-difference  techniques  allow  calculation 
in  two-  and  three-dimensional  structures.  In  general,  wave- theoretical 
techniques  applied  to  two-dimensional  synthetic  seismogram  calculations 
provide  highly  accurate  results  but  are  limited  in  their  practical  applica¬ 
tion  due  to  the  extensive  cost  and  computer  resources  required  for  their 
calculation.  In  contrast,  ray- theoretical  methods  provide  only  approxi¬ 
mate  calculations  and  are  limited  In  their  ability  to  simulate  wave  propaga¬ 
tion  effects,  but  are  computationally  efficient.  In  this  paper,  we  briefly 
describe  a  modification  of  the  disk-ray  theory  method  (Wiggins,  1976) 
for  calculation  of  synthetic  seismograms  in  two-dimensional  velocity 
structures  and  present  an  example  of  its  application  to  modeling  of  a 
complex  geologic  model  in  the  eastern  Snake  River  Plain  area  of  southern 
Idaho.  The  method  is  reasonably  accurate,  is  efficient  to  use  for  modeling 
purposes,  and  is  capable  of  calculating  seismograms  for  complex  two- 
dimensional  geologic  structures. 

A  variety  of  synthetic  seismogram  modeling  techniques  have  been 
developed  for  application  to  laterally  inhomogeneous  velocity  structures. 
Cerveny  et  al_.  (1977)  and  Cerveny  (1979)  present  a  theoretical  discussion 
and  results  of  model  studies  of  asymptotic  ray  theory  (ART)  calculations 
of  the  amplitudes  of  seismic  waves  propagating  through  an  inhomogeneous 
velocity  structure  approximating  the  waves  by  a  series  of  rays.  McMechan 
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and  Mooney  (1980)  described  a  modified  version  of  ART  and  showed  an  applica¬ 
tion  to  modeling  of  observed  crustal  seismic  data  for  the  Imperial  Valley, 
California.  Hong  and  Helmberger  (1978)  have  utilized  a  ray  theory  method 
known  as  glorified  optics  to  compute  synthetic  seismograms  for  primary 
and  multiple  reflections  for  some  homogeneous  layer  models  including 
non-planar  boundaries.  Frazer  and  Phinney  (1980)  have  described  the 
theoretical  background  for  a  generalized  phase  integral  method  which 
can  be  applied  to  generation  of  synthetic  seismograms  in  laterally  inhomo¬ 
geneous  media.  Wiggins  (1976)  and  Wiggins  and  Madrid  (1974)  have  developed 
a  ray  method  for  the  generation  of  synthetic  seismograms  in  which  the 
effects  of  wave  propagation  are  simulated  by  individual  rays  in  which 
the  propagation  can  be  thought  of  as  consisting  as  'disks'  of  energy 
traveling  along  the  raypath.  Each  disk,  which  is  perpendicular  to  the 
raypath,  intersects  the  surface  and  affects  an  area  surrounding  the  point 
of  emergence  of  the  ray.  Chapman  (1976)  described  a  first  motion  approxi¬ 
mation  which  is  formally  equivalent  to  the  disk  ray  theory  technique. 

All  of  the  above  mentioned  techniques  are  based  on  ray- theoretical 
solutions  to  the  elastic  wave  equation  in  two-dimensional  media.  They 
provide  only  approximate  solutions  for  particular  wave  propagation  phenomena 
and  are  limited  in  their  ability  to  simulate  wave  propagation  for  all 
types  of  wave  motion  and  under  all  geometrical  conditions  of  velocity 
structure.  For  example,  it  is  well  known  that  ray  theoretical  techniques 
are  inaccurate  in  the  area  of  a  caustic  or  turning  point  of  the  seismic 
rays  and  have  difficulty  in  the  generation  of  pure  head  waves.  Furthermore, 
ray  techniques* are  generally  limited  to  body  wave  propagation  problems. 

More  exact  wave  propagation  synthetic  seismogram  techniques  generally 
involve  a  direct  numerical  solution  to  the  elastic  wave  equation  in  two- 
dimensions.  Finite  difference  approximations  for  this  purpose  have  been 
used  by  Boore  (1972),  Alford  et  al_.  (1974)  and  Kelly  et  al_.  (1976) 
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and  finite  element  methods  are  described  by  Smith  (1975).  Because  these 
techniques  represent  complete  solutions  to  the  wave  equation,  the  finite 
difference  and  finite  element  synthetic  seismogram  methods  are  accurate 
and  provide  complete  seismograms  including  all  types  of  wave  phenomena 
in  both  homogeneous  and  inhomogeneous  velocity  structures.  However, 
large  requirements  of  computer  time  and  storage  for  these  techniques 
preclude  their  use  for  routine  modeling  of  two-dimensional  seismic  wave 
propagation. 

DISK  RAY  THEORY  METHOD 

The  disk  ray  theory  (DRT)  method  that  we  have  employed  follows  the 
theory  of  Wiggins  and  Madrid  (1974)  and  Wiggins  (1976).  They  presented 
a  theoretical  discussion  of  the  DRT  method  based  on  ray  approximations 
to  wave  propagation  in  laterally  homogeneous  media  and  pointed  out  the 
possibility  of  utilizing  such  techniques  for  laterally  inhomogeneous 
media  where  ray  tracing  results  would  provide  the  necessary  travel  time 
and  amplitude  information  for  constructing  DRT  synthetic  seismograms. 

The  DRT  method  that  we  have  utilized  is  briefly  described  below.  For 
a  discussion  of  the  theory  of  the  DRT  method  the  reader  is  referred  to 
Wiggins  and  Madrid  (1974)  and  Wiggins  (1976).  The  ray  tracing  data  required 
for  DRT  synthetic  seismogram  synthesis  follows  the  ART  techniques  of 
Cerveny  and  Ravindra  (1971)  and  Cerveny  et  aj_.  (1977). 

The  first  step  in  DRT  synthesis  for  two-dimensional  wave  propagation 
is  to  calculate  travel-time  and  amplitude  factors  for  seismic  raypaths 
propagating  for  phases  of  interest  through  the  velocity  structure.  For 
this  purpose,  we  utilize  the  'shooting  method'  of  ray  tracing  in  which 
rays  are  traced  by  iterative  application  of  the  generalized  Snell's  Law 
through  the  two-dimensional  velocity  structure  from  a  source  location 
given  an  initial  angle  of  incidence.  The  ray  may  refract  or  reflect 
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at  boundaries  and  may  Include  propagation  through  homogeneous  or  vertically 
and  laterally  inhomogeneous  layers.  The  velocity  structure  in  our  models 
is  given  by  layers  of  constant  or  slowly  varying  seismic  velocity  separated 
by  irregular  interfaces.  Alternatively  a  continuously  varying  velocity 
structure  with  vertical  and  lateral  gradients  can  be  approximated  by 
a  two-dimensional  grid  of  seismic  wave  velocities.  Attenuation  factors 
(Q~^)  are  also  given  for  the  model.  The  seismic  velocity  model  utilized 
employs  a  flat  earth  convention  and  thus  we  provide  for  small  positive 
velocity  gradients  within  homogeneous  layers  to  simulate  the  effects 
of  earth  curvature  according  to  the  earth  flattening  transformation  (Aki 
and  Richards,  1980,  p.  463-465).  Using  only  this  small  positive  velocity 
gradient  within  homogeneous  layers  corresponding  to  the  earth  flattening 
transformation  (apprbximately  0.001  km/s/km)  raypaths  corresponding  to 
pure  head  wave  propagation  can  be  generated. 

The  ray  tracing  through  the  two-dimensional  velocity  structure  provides 
travel  time,  distance,  raypath  parameter  (slowness)  and  an  amplitude 
factor  for  each  seismic  ray  for  input  into  the  DRT  synthesis  program. 

The  amplitude  factor  includes  the  amplitude  effects  of  geometrical  spreading, 
reflection  and  transmission  coefficients  at  interfaces,  a  free-surface 
conversion  factor  and  attenuation  due  to  absorption.  The  attenuation 

N  AT. 

is  calculated  from  an  exponential  amplitude  weighting  factor  (Exp(-Trf  z  --!)) 

i*l  4i 

according  to  the  travel-time  spent  within  each  Q  layer.  The  attenuation  due  to 
absorption  is  approximated  by  this  exponential  absorption  factor  corresonding 
to  the  dominant  frequency  of  the  wavelet  which  will  be  used  for  the  synthetic 
seismogram  generation.  Thus,  it  is  a  'single  frequency  Ql  approximation, 
and  results  in  attenuation  of  the  amplitude  of  propagated  waves,  but 
does  not  result  in  frequency  dependent  changes  in  waveform. 

A  flow  chart  shown  in  Figure  1  and  a  schematic  diagram  shown  in 
Figure  2  illustrate  the  process  of  ray  tracing  and  DRT  seismogram  synthesis 
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for  a  two-dimensional  velocity  model.  Ray  tracing  through  the  velocity 
structure  provides  time,  distance,  raypath  parameter  and  amplitude  factors 
(t,j,  x^ ,  p.  and  a.. )  which  can  be  used  to  generate  a  travel  time  curve 
for  all  of  the  phases  of  interest  which  are  included  in  the  ray  tracing. 

DRT  synthesis  sums  the  contributions  of  individual  amplitudes  for  raypaths 
in  the  vicinity  of  the  desired  seismogram  location  for  all  branches  of 
the  travel  time  curve  and  finally  convolves  the  amplitude  time  series 
with  a  desired  wavelet  to  produce  individual  seismograms.  Repeating 
the  amplitude  summation  for  all  distances  of  interest  produces  the  seismic 
record  section.  The  diagram  shown  in  Figure  2  illustrates  this  process 
for  a  single  seismogram  at  a  distance  Xq.  In  order  to  synthesize  this 
seismogram  by  the  DRT  method,  the  time,  distance,  amplitude  and  raypath 
parameter  values  for  the  individual  raypaths  are  generated  by  ray  tracing 
through  the  velocity  model.  For  each  seismic  phase  (branch  of  the  T-X 
curve),  the  amplitude  contributions  of  individual  raypaths  are  projected 
from  their  arrival  time  (using  the  slope  equal  to  the  raypath  parameter) 
to  the  distance  corresponding  to  the  seismogram  location  xQ.  Thus,  seismic 
energy  corresponding  to  raypaths  arriving  in  the  vicinity  of  the  distance 
Xq  contribute  to  the  amplitude  of  the  arrival  at  the  distance  Xq.  Distance 
and  time  weighting  factors  are  utilized  to  prevent  strong  arrivals  from 
large  distances  far  away  from  the  seismogram  location  being  projected 
into  the  seismogram  and  contributing  fictitious  amplitude.  Although 
the  distance  and  time  weighting  factors  are  chosen  arbitrarily,  we  have 
found  that  the  process  is  not  very  sensitive  to  the  choice  of  these  weighting 
factors.  The  ‘distance  factor  that  we  use  generally  corresponds'  to  about 
2  to  4  seismic  wavelengths  and  the  time  window  is  approximately  one  period 
of  the  dominant  frequency  of  the  wavelet.  After  the  individual  amplitude 
contributions  of  each  raypath  are  considered  for  each  of  the  possible 
branches  of  the  travel  time  curve  within  the  distance  given  by  the  distance 
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weighting  window,  the  individual  time  versus  amplitude  spikes  are  con¬ 
volved  with  the  wavelet  to  produce  the  DRT  synthetic  seismogram. 

MODEL  STUDIES 

In  order  to  test  the  DRT  method,  the  vertical  component  displacement 
theoretical  seismograms  (Figure  3)  were  computed  at  the  surface  of  a 
model  consisting  of  a  homogeneous  layer  over  a  half  space  for  an  explosive 
point  source  at  a  depth  close  to  the  surface.  The  result  is  compared 
with  modified  reflectivity  method  (MRM)  given  by  Kind  (1978). 

The  velocity  depth  model  was  chosen  arbitrarily  for  convenience 
so  that  all  the  major  arrivals  could  be  separated  clearly.  For  the  MRM 
seismograms,  shear  wave  velocities  and  densities  for  layers  of  the  model 
were  computed  from  the  "P-wave  velocities  assuming  a  Poisson's  ratio 
a  =  0.25.  Densities  were  computed  from  *  0.252  +  0.3788  (Birch, 

1964).  The  effects  of  attenuation  on  body  waves  is  included  by  introducing 
complex  velocities  (Braile,  1977). 

For  the  DRT  method,  the  Q"^  factor  which  is  dependent  on  a  single 
predominant  frequency,  surface  conversion  coefficients  (Cerveny  and  Ravindra, 
1971),  geometrical  spreading,  and  reflection  and  transmission  coefficients 
were  taken  into  account  in  the  amplitude  calculation.  For  simplicity, 
only  direct  waves  and  primary  P  wave  reflections  with  corresponding  head 
waves  are  considered  in  the  DRT  seismograms.  Later  arrivals  including 
multiple  reflections  are  identified  on  the  record  section  for  the  MRM 
seismograms.  Amplitudes  of  both  the  DRT  and  MRM  seismograms  were  multiplied 
by  distance  for  convenient  plot  scaling. 

Figure  3  shows  theoretical  seismograms  computed  by  both  the  MRM 
and  the  DRT  methods.  The  structure  is  a  homogeneous  layer  over  a  half 
space  with  only  the  earth  flattening  transformation  applied.  For  this 
case,  ray-theoretical  methods  generally  fails  to  predict  the  head  waves. 


but  in  our  DRT  computation,  the  head  wave  can  be  predicted  and  the  relative 
amplitude  and  waveform  character  agree  with  the  theoretical  seismograms 
computed  by  the  MRM.  The  major  difference  is  in  the  critical  point  region 
(near  10  km)  where  the  amplitude  for  the  DRT  method  is  larger  than  for 
the  MRM,  and  the  phase  is  distorted.  This  is  due  to  the  fact  that  wave 
theory  predicts  that  the  phase  change  due  to  reflection  varies  gradually 
in  the  near-critical  region,  whereas  ray  theory  predicts  that  the  phase 
shift  for  reflections  begin  abruptly  at  the  critical  distance.  Although 
minor  differences  are  apparent,  comparison  of  the  amplitude  and  waveform 
characteristics  for  the  DRT  and  MRM  synthetic  seismograms  shown  in  Figure 
3  indicates  that  the  DRT  method  is  capable  of  accurate  and  efficient 
synthetic  seismogram  calculations. 

APPLICATION  TO  MODELING  OBSERVED  SEISMIC  RECORD 
SECTIONS  FROM  THE  SNAKE  RIVER  PLAIN  AREA 

Sparlin  et  aK  (1982)  have  interpreted  a  136  km  reversed  seismic 
line  in  the  eastern  Snake  River  Plain,  Idaho  (SP  7  SE  in  the  northern 
Rocky  Mountains  province  to  the  Gay  Mine  shotpoint  located  in  the  Basin 
and  Range  province)  using  ray-trace  travel -time  modeling.  The  velocity 
structure  determined  by  Sparlin  et  al_.  (1982)  is  shown  in  Figure  4  along 
with  the  raypaths  of  principal  refracted  and  reflected  phases. 

In  this  section,  we  compute  the  corresponding  synthetic  seismograms 
for  the  Sparlin  et  al_.  (1982)  model  and  our  modified  model  using  DRT 
for  the  shotpoints  (SP  7  SE  and  Gay  Mine)  at  the  ends  of  the  velocity 
model.  Observed  and  synthetic  (DRT)  record  sections  for  the  Snake  River 
Plain  model  are  shown  in  Figures  5  and  6.  The  synthetic  seismograms 
calculated  for  the  model  proposed  by  Sparlin  et  a_K  (1982)  based  on  travel¬ 
time  modeling  (Figure  4)  display  amplitude  characteristics  for  phase  E 
(Figures  4,  5  and  6)  which  do  not  match  the  observed  data.  A  small 
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modification  to  the  Snake  River  Plain  velocity  model  (Figure  7)  results 
in  synthetic  seismograms  which  match  very  closely  the  amplitude  and  even 
waveform  character  of  the  observed  data  for  both  the  SP  7  SE  (Figure 
5)  and  Gay  Mine  (Figure  6)  shotpoints.  For  example,  the  A,  B,  C  phases 
(Figure  5)  show  the  expected  travel  time  delay  due  to  block  faulting. 

Large  amplitude  arrivals  are  calculated  near  65  km  in  the  D  phase  and 
from  60  to  110  km  for  the  F  phase  (Figure  5).  The  very  weak  amplitude 
of  phase  E  is  seen  on  the  synthetics  for  the  modified  model  as  observed 
on  the  real  data. 

In  order  to  illustrate  the  importance  of  two-dimensional  synthetic 
seismogram  modeling  such  as  for  the  Snake  River  Plain  example  (Figures 
5  and  6),  we  computed  MRM  synthetic  seismograms  for  a  plane-layered  model 
(laterally  homogeneous)  which  approximately  fits  the  travel  times  of 
th  observed  SP  7  SE  and  Gay  Mine  data.  Although  the  arrival  times  match 
reasonably  well,  the  amplitude  and  waveform  character  of  these  synthetics 
(Figure  8)  display  a  poor  comparison  to  the  observed  SP  7  SE  and  Gay 
Mine  record  sections.  The  excellent  match  of  the  two-dimensional  DRT 
synthetics  with  the  observed  data  and  the  inability  to  fit  the  observed 
data  with  a  laterally  homogeneous  model  provides  confirmation  of  the 
complex  velocity  model  beneath  the  eastern  Snake  River  Plain.  Additionally, 
this  example  illustrates  the  importance  of  amplitude  data  and  synthetic 
seismogram  modeling  for  detailed  and  accurate  interpretation  of  crustal 
velocity  structure. 


DISCUSSION  AND  CONCLUSIONS 

Comparison  of  the  synthetic  seismograms  computed  with  the  DRT  method 
with  those  of  the  reflectivity  method  have  shown  that  the  DRT  method 
is  efficient  and  reasonably  accurate.  The  method  can  be  applied  to  complex 
two-dimensional  velocity  and  Q  structures.  Application  to  modeling 
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observed  seismic  refraction  data  for  the  eastern  Snake  River  Plain  has 
illustrated  the  significance  of  two-dimensional  synthetic  seismogram 
techniques  and  the  utility  of  the  DRT  method.  Modern  refraction /wide- 
angle  reflection  data  for  crustal  seismic  studies  will  necessitate  increased 
use  of  modeling  techniques  capable  of  laterally  inhomogeneous  velocity 
structure.  Although  the  DRT  method  is  an  approximate  technique  and  is 
some  what  restricted  in  its  application,  its  efficiency  and  accuracy 
make  it  suitable  for  routine  modeling  of  data  in  laterally  inhomogeneous 
applications.  The  principal  limitation  that  we  have  found  with  the  DRT 
technique  is  the  inability  to  trace  rays  through  some  complex  geologic 
structures.  This  restriction  is  more  a  limitation  of  the  ray  tracing 
method  than  of  the  DRT  seismogram  synthesis. 
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FIGURE  CAPTIONS 


Figure  1.  Flow  chart  illustrating  the  process  of  calculation  of  synthetic 
seismograms  by  the  disk-ray  theory  method. 

Figure  2.  Schematic  diagram  illustrating  the  calculation  of  synthetic 

seismograms  by  the  disk-ray  theory  method.  Upper  diagram  shows 
a  travel  time  curve  with  a  seismogram  synthesized  at  the  distance 
xq.  Middle  diagram  illustrates  the  velocity  model  with  compres- 
sional  wave  velocities  (a)  and  attenuation  factors  (Q)  described 
by  a  two-dimensional  model  including  non-planar  interfaces. 

Raypaths  for  direct,  reflected  and  refracted  waves  through 
the  velocity  model  are  illustrated.  Lower  diagram  illustrates 
the  process  of  DRT  amplitude  synthesis  in  which  the  amplitudes 
for  rays  for  a  given  travel  time  branch  are  added  by  projecting 
their  time  of  arrival  along  the  travel  time  curve  and  applying 
both  distance  and  time  weighting  factors  to  the  amplitude  con¬ 
tributions  of  each  raypath.  The  response  at  the  distance  xQ 
is  the  sum  of  the  weighted  responses  of  each  raypath  for  each 
travel -time  branch  projected  to  the  appropriated  arrival  time. 
Convolution  of  the  individual  projected  ray  amplitudes  with 
a  source  wavelet  completes  the  disk  ray  theory  seismogram  synthesis. 

Figure  3.  Synthetic  seismogram  record  section  calculated  using  the  modified 
reflectivity  method  (MRM)  and  disk  ray  theory  (DRT)  method  for 
a  plane  homogeneous  layer  overlying  a  half-space.  The  single 
layer  is  5  km  thick  with  a  velocity  of  4.5  km/s  overlying  ?. 
half-space  of  velocity  6.8  kra/s.  Q  values  are  150  and  500  for 
the  layer  and  half-space  respectively.  The  DRT  seismograms 
are  displaced  slightly  to  the  left  of  the  correct  distance  in 
the  plot  for  clarity.  The  ray  diagrams  at  the  bottom  of  the 
figure  illustrates  the  notation  for  the  travel  time  branches. 

Figure  4.  Crustal  model  interpreted  by  Sparlin  et  al_.  (1982)  for  the  eastern 
Snake  River  Plain  in  Idaho.  The  raypaths  for  a  variety  of  refracted 
(upper  diagram)  and  reflected  (lower  diagram)  phases  are  shown. 
Numbers  are  compressional  wave  velocities  in  km/s.  The  datum 
is  1.2  km  elevation.  The  location  of  shotpoint  7  and  the  Gay 
Mine  shotpoint  are  at  the  northwest  and  southeast  ends  of  the 
model,  respectively.  NRM  is  northern  Rocky  Mountains  province; 

SRP  is  the  Snake  River  Plain  province;  and  BR  is  the  Basin  and 
Range  province. 

Figure  5.  Observed  and  synthetic  seismograms  for  the  SP  7  SE  profile  across 
the  velocity  structure  shown  in  Figure  4.  Upper  diagram  shows 
the  synthetic  seismograms  calculated  by  the  DRT  method  for  the 
Sparlin- et  ajL  (1982)  model  (Figure  4).  Middle  diagram  illus¬ 
trates  tTTe  amplitude-corrected  observed  data  for  the  SP  7  SE 
profile  recorded  during  the  Y-SRP  1978  experiment.  Lower  diagram 
shows  the  synthetic  seismograms  calculated  by  the  DRT  method 
for  the  modified  SP  7  SE  model  (Figure  7).  The  amplitudes  for 
all  of  the  seismograms  shown  in  this  figure  are  scaled  with 
a  distance  factor  of  amplitudes  times  distance  to  the  1.5  power 
for  convenient  amplitude  scaling.  The  phase  notation  shown 
In  the  middle  seismic  section  corresponds  to  the  phases  illus¬ 
trated  in  Figure  4. 


72 


Figure  6.  Observed  and  synthetic  seismograms  for  the  Gay  Mine  shotpoint 
(Figure  4).  Upper  diagram  illustrates  the  synthetic  seismo¬ 
grams  calculated  by  the  DRT  method  for  the  Gay  Mine  shotpoint 
and  the  velocity  structure  shown  in  Figure  4.  Middle  diagram 
shows  observed  data  from  the  Gay  Mine  shotpoint.  Lower  diagram 
illustrates  synthetic  seismograms  calculated  by  the  DRT  method 
for  the  modified  Gay  Mine  model  (Figure  7).  All  of  the  seismo¬ 
grams  in  these  record  sections  have  been  scaled  with  the  identical 
scaling  factor  of  amplitude  times  distance  to  the  1.5  power 
to  provide  for  convenient  plot  scaling. 

Figure  7.  Velocity  structure  for  the  modified  eastern  Snake  River  Plain 
model.  The  principal  modification  of  this  model  as  compared 
to  the  Sparlin  et  al_.  (1982)  model  illustrated  in  Figure  4  is 
a  flattening  of  the  Interface  separating  the  6.53  km/s  layer 
from  the  6.15  km/s  layer  beneath  the  Snake  River  Plain  (SRP). 

This  modification  to  the  velocity  structure  primarily  affects 
calculation  of  the  phase  E  (Figure  4)  for  which  raypaths  are 
illustrated.  Numbers  in  parentheses  are  Q  values  used  in  the 
DRT  calculations  for  the  SP  7  SE  and  Gay  Mine  synthetic  seismo¬ 
grams  (Figures  5  and  6)  and  are  based  on  the  Q  values  interpreted 
by  Braile  et  a^.  (1982). 

Figure  8.  Synthetic  seismograms  calculated  by  the  modified  reflectivity 

method  for  a  plane  layered  structure  whose  travel  times  approxi¬ 
mate  the  observed  travel  time  data  for  the  SP  7  SE  and  Gay  Mine 
observed  record  sections.  The  velocity  and  Q  structure  are 
illustrated  In  the  left  hand  diagram.  Amplitudes  of  the  seismo¬ 
grams  are  multiplied  by  distance  to  the  1.5  power  for  convenient 
plot  scaling. 
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